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Battery  cell  life  depends  critically  on  how  the  cell  is  used.  Therefore,  battery  chargers  and  battery  manage¬ 
ment  systems  must  be  designed  to  control  cell  usage  carefully.  In  order  to  design  optimal  battery  controls 
that  effect  a  tradeoff  between  cell  performance  (in  some  sense)  and  cell  life,  a  model  of  cell  degradation 
is  necessary.  This  model  must  be  simple  and  incremental  in  order  to  be  implemented  by  an  inexpensive 
microcontroller.  This  paper  takes  a  first  step  toward  developing  such  a  controls-oriented  comprehensive 
cell  degradation  model  by  deriving  a  reduced-order  model  of  a  single  mechanism:  lithium  deposition  on 
overcharge,  along  with  the  resulting  resistance  rise  and  capacity  loss.  This  reduced-order  model  approxi¬ 
mates  a  physics-based  PDE  model  from  the  literature,  is  simple  and  accurate,  and  may  be  used  in  optimal 
strategies  for  controlling  lithium  ion  batteries. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Battery  packs  require  monitoring  and  control  to  ensure  safety, 
to  deliver  peak  performance,  and  to  maximize  life.  It  is  gener¬ 
ally  necessary  to  make  compromises  among  these  objectives  since 
high  charge  and  discharge  rates  accelerate  aging  and  can  result 
in  unsafe  operating  conditions.  Accordingly,  a  present  goal  of  our 
research  program  is  to  design  battery  control  methodologies  to 
mathematically  optimize  a  tradeoff  between  performance  and  life. 
To  do  so,  we  require  quantitative  models  of  the  dominant  cell 
degradation  mechanisms  that  can  accurately  predict  a  measure 
of  the  degradation  that  would  be  caused  by  candidate  control 
actions.  These  models  must  be  simple  enough  to  be  executed 
quickly  on  an  inexpensive  embedded  systems  processor,  so  we 
seek  to  develop  reduced-order  approximations  to  more  complex 
partial-differential  equation  (PDE)  models  describing  mechanical 
and  chemical  cell  degradation  processes. 

Presently,  many  battery  management  systems  track  macro¬ 
scopic  indicators  of  aging  such  as  capacity  fade  and  power  fade  but 
few,  if  any,  track  physical  indicators  of  aging  such  as  degree  of  solid- 
electrolyte-interphase  (SEI)  layer  growth  on  anode  particles,  degree 
of  lithium  plating  on  anode  particles,  degree  of  cathode  dissolution, 
and  so  forth.  However,  identical  combinations  of  capacity  loss  and 
resistance  rise  might  be  achieved  by  traveling  different  paths— by 
exciting  different  physical  mechanisms  of  degradation— yielding 
potentially  different  resultant  optimal  control  strategies  for 
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operating  the  battery  pack.  Therefore,  knowing  cell  resistance  and 
capacity  alone  is  not  sufficient  to  extend  cell  life  to  the  maximum 
extent  via  control. 

We  believe  that  by  modeling  the  dominant  physical  degradation 
mechanisms  that  occur  in  a  cell,  and  then  by  using  those  mod¬ 
els  in  an  optimized  predictive  control  algorithm,  battery  life  can 
be  extended.  For  a  practical  implementation,  this  requires  that  the 
degradation  model  be  simple  and  incremental.  We  consider  a  model 
to  be  simple  if  it  is  described  by  a  small  number  of  ordinary  differ¬ 
ence  equations,  and  incremental  (or  recursive)  if  it  predicts  a  future 
degradation  state  based  on  a  present  degradation  state  and  a  pro¬ 
posed  cell  current.  The  model  will  necessarily  be  a  function  of  cell 
temperature,  state-of-charge  (SOC),  and  possibly  other  factors  as 
well. 

Degradation  leading  to  capacity  loss  and  resistance  rise  can 
occur  either  due  to  mechanical  stress  factors  or  chemical  side 
reactions  [1-8].  Various  schemes  have  been  introduced  to  model 
side  reactions,  including  that  of  Darling  and  Newman,  who  first 
introduced  the  concept  of  modeling  parasitic  effects  in  lithium-ion 
batteries  by  modeling  solvent  oxidation  side  reactions  [3].  Here,  we 
focus  on  overcharge,  and  on  creating  a  simplified  approximation  to 
a  model  proposed  by  Arora  and  colleagues  [9].  Overcharge  mani¬ 
fests  first  as  a  metallic  lithium  deposit  on  the  surface  of  the  negative 
electrode  solid  particles  during  charge,  predominantly  near  the 
separator.  Subsequently,  the  lithium  can  and  does  further  combine 
with  electrolyte  material  to  form  other  compounds  such  as  Li20, 
LiF,  Li2C03,  and  polyolefins.  The  nature  of  the  final  product  is  not 
our  major  concern  in  this  paper;  rather,  the  issue  is  that  lithium  is 
irreversibly  lost.  This  phenomenon  is  an  unintended  side  reaction 
that  can  lead  to  severe  capacity  fade,  electrolyte  degradation,  and 
a  possible  safety  hazard. 
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The  present  paper  builds  on  and  extends  the  work  in  [9]  to  pro¬ 
pose  a  simple  incremental  model  of  overcharge  and  associated 
capacity  loss  and  resistance  rise.  The  proposed  model  may  then 
be  used  in  an  optimal  control  scheme  that  will  reduce  side  reac¬ 
tions  and  limit  cell  degradation.  (The  control  task  is  the  subject  of 
present  research,  and  will  be  reported  in  a  later  paper.  The  method¬ 
ology  presented  in  [10]  is  one  possible  approach.)  While  various 
efforts  have  been  made  to  create  reduced-order  models  of  ideal-cell 
dynamics  (i.e.,  models  that  don’t  consider  degradation),  includ¬ 
ing  single  particle  models  by  Chaturvedi  [11]  and  one-dimensional 
models  by  Smith  [12],  we  believe  this  to  be  one  of  the  first  attempts 
to  create  a  reduced-order  model  to  describe  a  side  reaction.  This 
reduced-order  model  of  overcharge  is  a  first  step  toward  creating  a 
complete  coupled  reduced-order  model  of  all  dominant  cell  degra¬ 
dation  mechanisms,  which  then  could  be  used  in  an  optimal  control 
scheme. 

This  paper  will  proceed  by  reviewing  the  model  of  overcharge 
proposed  in  [9].  A  reduced-order  one-dimensional  approximate 
model  is  then  derived.  Results  of  pulsed-current  simulations  of  the 
full-order  model  and  the  reduced-order  model  are  presented  and 
compared,  results  are  discussed,  and  conclusions  are  made. 


2.  Physics-based  model  of  overcharge 

Changes  at  the  electrode/electrolyte  interface  due  to  side  reac¬ 
tions  between  the  anode  and  electrolyte  are  considered  to  be 
one  of  the  primary  causes  of  anode  aging  [8].  One  particularly 
severe  degradation  mechanism  is  that  of  overcharge,  where  metal¬ 
lic  lithium  plates  on  the  electrode  surface.  A  physics-based  coupled 
PDE  model  of  ideal-cell  and  overcharge  dynamics  has  been  reported 
in  [9].  Our  goal  is  to  create  a  high-fidelity  reduced-order  model  of 
this  PDE  degradation  model;  therefore,  we  adopt  the  same  assump¬ 
tions  as  they,  which  were: 


1.  The  main  side  reaction  is  expressed  as  Li+  +  e_  ->  Li(s),  which 
occurs  at  Us  =  0  V  versus  Li/Li+  during  an  overcharge  event.  This 
lithium  metal  is  expected  to  form  first  near  the  electrode¬ 
separator  boundary  where  the  surface  overpotential  is  greatest. 

2.  Lithium  metal  deposited  on  the  negative  electrode  reacts  quickly 
with  solvent  or  salt  molecules  in  the  vicinity,  yielding  Li2C03,  LiF, 
or  other  insoluble  products.  A  thin  film  of  these  products  protects 
the  solid  lithium  from  reacting  with  the  electrolyte.  Solid  lithium 
can  still  dissolve  during  discharge,  but  once  lithium  is  consumed 
in  a  insoluble  product,  it  is  permanently  lost. 

3.  The  insoluble  products  formed  are  a  mixture  of  different  species, 
resulting  in  averaged  mass  and  density  constants  used  in  a  later 
equation  describing  the  formation  and  growth  of  a  resistive  film. 

4.  Only  the  overcharge  reaction  is  considered  (e.g.,  SEI  film  growth 
and  other  degradation  mechanisms  are  not  modeled). 


The  overcharge  model  of  [9]  is  tightly  coupled  with  a  “pseudo 
two-dimensional”  Doyle-Fuller-Newman  porous-electrode  style 
model  of  ideal-cell  dynamics  [4],  which  assumes  that  the  solid  and 
electrolyte  phases  are  considered  continuous  and  gives  no  consid¬ 
eration  to  the  underlying  microstructure  of  the  electrode.  While 
porous-electrode  theory  is  well  known,  there  are  variants  in  its 
implementation,  so  we  present  the  full  set  of  model  equations 
here,  in  summary  form.  We  list  the  model  in  terms  of  the  negative 
electrode  only;  all  equations  have  exact  analogs  for  the  positive 
electrode  as  well,  with  the  exception  of  the  lithium  plating  equa¬ 
tions,  as  plating  happens  only  on  the  negative  electrode.  Symbols 
are  defined  in  a  table  at  the  end  of  this  paper;  only  the  most  relevant 
are  identified  here. 


The  conservation  of  lithium  in  a  single  particle  is  modeled  by 
Fields  law  of  diffusion 

3 cs(r,  x,  t)  Ds  3  /  2  9cs(r,  x,  t) 

3 1  r2  dr  V  3 r 


with  boundary  conditions 
3 cs(r,  x,  t) 


dr 


=  0  and  D« 


3 cs(r,  x,  t) 


dr 


-jn(x,  t) 
asF  ’ 


where  jn(x,  t)  is  the  intercalation  current  density.  Conservation  of 
lithium  in  the  electrolyte  phase  gives  the  following  equation, 


3ce(x,  t)  _  3  /  Dfff  3ce(x,  t)\  (1  -  )  die(x,  t) 

dt  dx  \  Se  dx  )  +  Fse  dx 


(1) 


where  Dff  =  Des^rug  and  having  the  following  zero  flux  boundary 
conditions  at  the  current  collectors, 


3ce(x,  t) 


dx 


3 ce(x,  t) 


x=0 


3x 


=  0. 


x=L 


In  Eq.  (1 ),  the  current  density  in  the  electrolyte  has  the  form, 
3 ie(x,  t) 


dx 


—  j total (X,  t), 


(2) 


(3) 


where  jtotai(x,  t)  =jn(x,  t)  +js(x,  t)  andjs(x,  t)  is  the  side-reaction  cur¬ 
rent  density.  Eq.  (3)  has  the  boundary  conditions 


ze(0,  t)  =  0  and  ie{Ln,t)  = 


m 

A  ’ 


(4) 


where  J(t)  is  the  cell  current  and  A  is  the  current-collector  area. 
Current  density  in  the  solid  has  the  form 


is(x,  t)  =  ^  -  ie(x,  t ). 


(5) 


Charge  conservation  in  the  solid  phase  in  the  electrodes  is 
described  by  Ohm’s  law, 


-<7eff  0s(x,  t)  =  is(x,  t), 
dx 


(6) 


where  creff  =  ers^' rug,  <ps{ 0,  t)  =  0  by  convention,  and  having  boundary 
conditions  at  the  current  collector  and  separator  boundaries, 


eff  d<t>s{x,  t ) 


dx 


and 

A 


3 0s(x,  t) 


dx 


=  0. 


x=Ln 


Conservation  of  charge  in  the  electrolyte  gives  the  equation 

3 0e(x,  t)  -ie(x,  t)  x^ff  3  In  ce(x,  t) 


dx 


-eff 


-eff 


dx 


(7) 


(8) 


where  /ceff  =  Ksfug,  /c^ff  =  (2 RgT/ceff/E)(t°  -  1),  and  having  bound¬ 
ary  conditions 


30e(*,  t) 


dx 


30e(X,  t) 


x=0 


dx 


=  0. 


x=L 


The  PDEs  are  coupled  via  the  intercalation  current  density, 
expressed  as  the  Butler-Volmer  equation, 


jn(x,  t)  =  a„i0,„ 


exp 


Ota.nF 

~w' 


r]n{x,  t)  -  exp  - 


Otc.nF 

w 


r]nix,  t ) 


(9) 


which  is  driven  by  the  overpotential 
Ux,  t )  =  4>s(X,  t )  -  0e(x,  t)-Un- 
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where  R®El  is  the  initial  film  resistance  that  is  produced  during  the 
formation  period  of  the  battery. 

In  addition  to  the  resistance  change,  there  is  a  capacity  loss 
caused  by  the  side  reaction  current  during  charge,  leading  to  capac¬ 
ity  changing  via  the  relationship 

^l=A^njs{x,t)  dx.  (13) 

3.  Derivation  of  reduced  order  model  and  plating 
prediction  equations 


where  i0,n  is  the  exchange  current  density, 

<0,n  =  /<n(csmnax  -  Cs,n)aa,n {Cs,n)ac,n {ce)°la,n , 

and  Un  is  the  equilibrium  potential  which  is  evaluated  as  a  function 
of  the  solid  phase  concentration  at  the  surface  of  the  particle. 

With  this  “ideal-cell”  (no  degradation)  model  in  place,  it  is 
straightforward  to  add  the  terms  that  account  for  lithium  plat¬ 
ing.  The  side-reaction  current  density  js  (i.e.,  the  rate  of  irreversible 
lithium  loss  due  to  lithium  plating)  can  be  expressed  as 


To  effect  an  optimal  control  strategy,  the  processors  in  battery 
chargers  and  battery  management  systems  must  be  able  to  calcu¬ 
late  the  rate  of  lithium  loss  due  to  overcharge  very  quickly  and 
accurately.  Solving  the  coupled  PDE  equations  described  above  is 
too  complicated  for  such  a  process.  The  degradation  model  needs 
to  be  much  faster  and  simpler.  In  this  section,  we  present  a  simpler 
incremental  model  for  calculating  js(x,  t),  Rr lm(t),  and  Q(t),  which 
yields  an  algebraic  rather  than  differential  solution. 

To  create  the  reduced-order  model,  several  additional  assump¬ 
tions  are  made: 


js(x,  t)  =  min  (V  ani0,s  exp  ^ 


RST 


W, 


«c,sf 


RgT 


-exp 

where  aa,s  =f=  oic,s  in  general 


*7s(x,  t ) 


(10) 


T]s{X,  t )  =  0S(X,  t)  -  0e(X,  t)  —  Us  — 


js{X ,  t) 

On 


and  where  the  side-reaction  exchange  current  density 
io,s  =  1<n,s{Ce)aci’s .  The  side  reaction  is  semi-irreversible  in  the 
sense  that  it  includes  an  anodic  rate  term,  but  does  not  allow  an 
overall  positive  side-reaction  current  density.  The  side  reaction 
occurs  only  at  spatial  locations  in  the  anode  where  r]s{x ,  t)  <  0.  This 
is  enforced  in  Eq.(10)  by  the  “min”  function,  which  sets  js(x,  t)  =  0 
for  values  of  x  where  r]s{x,  t)  >  0,  but  to  the  value  computed  by 
the  Butler-Volmer  equation  when  r]s{x ,  t)<0.  A  typical  scenario 
is  plotted  in  Fig.  1,  where  rjs{x,  t )  is  sketched  across  the  electrode 
width.  In  this  example,  plating  will  occur  in  the  interval  from  x  =  x0 
to  x  =  Ln.  Note  that  this  illustration  shows  that  the  cell  can  be  quite 
far  away  from  100%  state-of-charge  and  still  have  plating  occur 
near  the  separator  if  a  large  enough  charge  current  pulse  is  applied 
to  the  cell’s  terminals.  The  state-of-charge  is  only  one  variable 
of  importance— ultimately,  the  local  overpotential  determines 
whether  plating  occurs. 

Once  the  local  side  reaction  current  has  been  calculated,  the 
change  in  the  film  thickness  8f[\m  during  charging  is  modeled  as 


^film(x?  t) 

3 t 


M 

anpF 


js{x,  t), 


(11) 


where  M  and  p  are  the  average  molecular  weight  and  density  of 
lithium  and  products.  The  resistance  ^products  (x. 0  of  the  film  formed 
by  lithium  and  other  products  (Li2C03  is  used  as  an  example  here) 
is  given  by 


^products (*’  0-Zli  +ZU2C03 

where  zLi  and  zLi2C03  are  the  volume  fractions  of  lithium  and  Li2C03 
present  in  the  film.  The  film  resistance  is  then  given  by 

<^film(x?  0  —  ^SEI^’  ^)  ^products(X?  0  (12) 


1.  The  cell  is  always  in  a  quasi-equilibrium  state,  allowing  the 
exchange  current  density  z0,n  to  be  calculated  from  the  cell  SOC 
alone,  neglecting  local  variations  in  electrolyte  and  solid  surface 
concentration.  The  estimated  value  of  js(x,  t)  then  corresponds 
to  a  suddenly  applied  current  pulse  of  magnitude  J(t),  which  is 
constant  over  some  time  interval  At. 

2.  The  total  current  density  is  uniform  over  the  anode:  jtotai(x,  t)  = 
./total (0-  This  assumption  decides  the  solution  profile  of  other 
variables.  The  solution  profile  inside  of  batteries  will  be  more 
complicated,  and  the  results  section  investigates  the  accuracy 
lost  by  making  this  assumption. 

3.  We  model  the  short-term  time  rate  of  change  of  concentration  of 
lithium  in  the  electrolyte  as  proportional  to  the  applied  current 
spread  out  over  the  electrode:  3ce(x,  t)/3t« / 3I{t)l{FALn ). 

4.  We  assume  that  the  absolute  variation  in  lithium  concentration 
in  the  electrolyte  is  negligible  in  comparison  to  its  gradient.  This 
allows  us  to  write  3lnce(x,  t)/3x  =  (l/ce(x,  t))(3ce(x,  t)/3x)  « 
(l/ce)(3ce(x,  t)/3x),  where  ce  is  the  volume-average  concentra¬ 
tion  of  lithium  in  the  electrolyte. 

5.  The  anodic  and  cathodic  charge-transfer  coefficients  of  the  inter¬ 
calation  kinetics  are  equal  (an,a  =  &n,c  =  0.5). 


With  these  models  and  assumptions  in  mind,  we  take  the  fol¬ 
lowing  steps  to  derive  the  rate  of  lithium  loss  due  to  overcharge: 


1.  We  derive  the  ionic  current  ze(x,  t)  by  integrating  Eq.  (3). 

2.  We  show  that  0s(x,  t) «  0  at  all  points  in  the  negative  electrode. 

3.  We  substitute  ie(x,  t)  into  Eq.(l)  to  determine  the  electrolyte 
concentration  gradient  3ce(x,  t)/3x. 

4.  We  substitute  ze(x,  t )  and  3ce(x,  t)/3t  into  Eq.  (8)  to  determine  the 
electrolyte  potential  0e(x,  t). 

5.  Finally,  we  determine  the  rate  of  lithium  loss  due  to  overcharge 
from  0s(x,  t)  -  0e(x,  t). 


3.1.  Determining  the  ionic  current 

To  determine  the  ionic  current,  we  integrate  Eq.  (3)  across  a  por¬ 
tion  of  the  negative  electrode.  We  employ  Liebnitz’  integration  rule 
to  correctly  deal  with  the  variable  integration  limits,  resulting  in: 
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r^dx  -fwx.t  >dx 

Jo  9X  Jo 

ie(x,t)~  ie(0^t)  =jtotal(t)x 
0,  by  Eq.(4) 

te(X,t)  =] total(f)X,  0  <X  <Ln. 

Since  ie(Tn,  t)  =  J(t)M  by  Eq.(4),  we  have  that  AotaiW  =I(t)/(ALn). 
With  this  substitution,  we  arrive  at 

ie(x,  t )  =  £x.  (14) 

Note  that  the  sign  of  /(t)  is  negative  during  charge,  and  positive 
during  discharge. 

3.2.  Showing  that  4>s(x,  t)^0 


90e(*,  t) 
dx 


-ie{x,t)  d In ce(x,  t) _  -/(t)  /c^ff  d ce(x,  t) 
/ceff  /ceff  dx  /ceffALn  X  /ceffce  dx 


4ff^e-(l-t+)  1  /(t) 

Ce  D|ffF  J  Keff/4In  ’ 


0  <x  <Ln. 

(18) 


The  gradient  of  the  electrolyte  potential  is  linear  in  x,  meaning  that 
the  electrolyte  potential  itself  is  parabolic  in  x.  This  is  again  a  very 
good  approximation  to  the  full  order  model  of  a  cell. 

We  will  integrate  once  again  to  determine  </>e(x,  t).  However,  we 
first  clean  up  notation  by  defining 


E(t)  = 


Kjf  fee  -  (1  -  fp  )  1  Ijt) 

Ce  DfF  \  KeffALn ' 


(19) 


Then,  d<pe{x,  t)/9x  =  E(t)x.  We  integrate  this  expression  to  find  <pe(x, 
t) 


We  can  now  solve  for  the  solid  potential.  Integrating  Eq.  (6)  and 
employing  Eqs.  (5)  and  (14)  gives 


/ 


~a  g^0s(X,t)d/ 

<ps{X,  t )  -  4>s[ 0,  t) 


is(x>  Od/ 


j_  r  Kt) 
Wo  A 


4>six ,  t) 


M.  \X_*L 

aeffA  2  ln 


0  <x<Ln.  (is) 


The  value  of  0s(x,  t )  is  later  used  in  determining  whether  plating 
occurs  by  comparison  with  (j)e(x,  t).  Due  to  the  high  conductivity  of 
the  solid  versus  the  electrolyte,  we  find  that  4>s{x )  is  in  the  order  of 
microvolts,  while  <fie{x,  t)  is  in  the  order  of  millivolts.  Therefore,  we 
proceed  by  assuming  0s(x,  t)^0.  However,  Eq.(15)  could  be  used 
instead. 


3.3.  Determining  the  electrolyte  concentration  gradient 


rmx>dx  ,rmx  dx 

Jo  Jo 

(peix,  t)  -  0e(O,  t)  =  ^X2 
E(t) 

0e(X,  t)  =  -^X2  +  0e(O,  t),  0  <  X  <  Ln.  (20) 

To  complete  the  calculation  of  </>e(x,  t),  we  must  determine  </>e(0, 
t).  To  do  so,  we  select  the  value  of  0e( 0,  t)  tha^will  result  in  the 
average  0e(x,  t )  across  the  length  of  the  anode  </>e(t)  being  equal  to 
the  expected  average.  The  expected  average  is  found  by  inverting 
the  Butler-Volmer  Eq.  (9),  and  by  finding  the  value  of  77  that  results 
in jn.  This  is  easily  done  if  aa,n  =  oic,n  =  0.5  (cf.  Assumptions),  when 

*?(t) =  — jf-asinh  f  .  (21) 

F  \2anio,n  J 

Then,  fe{t)  =  -(77(f)  +  Un  +]n{t)Rmm/an). 

To  use  this  result  to  determine  0e(  0,  t),  we  must  calculate  the 
corresponding  value  of  0e(t)  from  Eq.  (20). 


It  is  necessary  to  correlate  the  applied  current  to  the  change  of 
concentration  in  the  electrolyte.  Integration  of  the  concentration 
dynamics  Eq.  ( 1 )  across  a  portion  of  the  negative  electrode,  applying 
Assumption  3,  gives 


/ 


dce(X ,  t ) 
dt 


dcj(x,  0 
dx2 


die(X,  t) 

dx 


dx 


3ce(x,  t) 
3x 


yg/(f)  _  Df 
FALn  X  ~  Se 


dce(x,  t )  3ce(0,  t) 


dx 


foe -(!-#)/(<■) 
DfF  AL„ 


1  n-ft)j(o 

Fee  ALn 


x,  0  <  x  <  Ln. 


(16) 

(17) 


Note  that  the  gradient  of  concentration  is  linear  in  x,  thus  implying 
the  concentration  itself  is  parabolic  in  x  across  the  electrode.  We 
find  this  to  be  a  very  good  approximation  to  the  continuum  model 
of  a  cell. 


3.4.  Determining  the  electrolyte  potential 


Mt)  = 


if 


0e(x,  0  d*  =  p~ 


—J—X^  +  0e(O,  t)x 

0 


’  =  ^n+0e(O,  t ) 


Pei 0,  t )  =  0e(t)  -  =  -  (  ^-asinh 


2RrT 


Jn(0 

2anio  n 


-ll. 


Finally,  this  gives  us  the  equation 


Fit).  ( 2RgT 


Peix ,  t)  =  ^X2  - 

_m,  2 

6 


asinh 


Jn(f) 

2dnio,n 


+  un  + 


+  Un  + 


jn(t)R  filn 


(22) 


Jn(f)^filn 


3nce  again,  to  ease  notation  we  define 

Wasinh(ie 

2anlQ  n 


F 


+  Un  + 


Jn  (O^filn 

an 


which  allows  us  to  write 
Mx,  r)Wx2-P(t). 


(23) 


(24) 


With  concentration  of  the  electrolyte  derived,  development  of 
the  electrolyte  potential  can  proceed.  Substitute  Eqs.  (17)  and  (14) 
into  Eq.  (8)  and  employ  Assumption  4  to  derive  the  potential  of  the 
electrolyte  across  the  electrode: 


3.5.  Plating  equations 

Lithium  plating  occurs  only  during  charge,  and  at  any  point 
along  the  anode  where  77 s(x,  t)<0  (by  Eq.(10)).  As  we  are 
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assuming  that  0s(x,  t)«  0  and  Us  =  0,  plating  occurs  where  0e(x, 
t)  +  (js(x,  t)Rfi\mlan )  >  0.  The  parabolic  shape  of  Eq.  (24)  means  that  if 
0e(O,  t)  >  —js(x,  t)Rfilm/an  during  charge,  there  is  plating  across  the 
entire  width  of  the  anode.  If  0e(Ln,  t)  <  -  js(x ,  t)Rfiim/an,  then  there  is 
no  plating  at  any  location  across  the  width  of  the  anode.  Otherwise, 
there  is  plating  that  occurs  across  a  region  from  x0  to  Ln,  where  x0 
is  a  value  yet  to  be  determined. 

To  find  the  location  where  plating  begins,  we  solve  for  0e(*o,  0  = 


5.  Compute 


0, 


x0  =  <  Ln, 


,P(t)  Js(f)^film/^n 


E(t) 


<  JsWR film  . 

On 

E(t)l2  <  2 P(t)  -  2Js(f)Rfilm 

an 

otherwise. 


-j-Xq  =P(t)- 


x0  =  \  2 


On 


,r(t)  Js(f)Pfilm/^n 


E(t) 


6.  Compute 


*7s,oc(0—  t 
Ln 


L3n  -xg)-P(t)(I„-x o)+^(^Rfiln 
o  an 


(Ln  —  Xq) 


We  “spread  out”  the  plating  region  over  the  electrode  width  to  get 
an  equivalent  overpotential, 


WO  =  T~  /  — 0e(x,  t)  -  dx 


1  fU 

L»  k  ' 


J2  Df  I  7s(0^filn 


=  --  /  ^  -P(0  + 

Jx0 


On 


dx 


1 

Tn 

1 

Tn 


E[px3-P(t)x  +  js(t)Rfi"1 
b  an 


*0 


v3 ^  ^  ,  JsCO^filn 


^(^-x^)-P(t)(Ln-Xo)+. 
o  an 


(Ln  —  Xo)  • 
(25) 


The  side-reaction  deposition  current  from  overcharge  is  also 
described  by  the  Butler-Volmer  equation,  except  in  this  case  the 
anodic  and  cathodic  transfer  coefficient  are  typically  not  equal  and 
the  reaction  is  favored  to  be  cathodic: 


Js(0  =  ank,s 


exp 


f  Oia,sF 

\Rst 


-  exp 


(  oic,sF 

V  RST 


(26) 


7.  Compute 


m  = 


1  anio,s 

_exp  (^ ) 

-exp  (-if ) 

lo, 

Xo  <  Ln 
Xo  =L„. 


8.  Iterate  steps  3-7  until  \jn(t)  +  js(t)  -  ( I(t)/(ALn))\  <  €,  where  6  is 
a  convergence  factor. 

The  iteration  in  step  8  is  generally  done  using  a  nonlinear  solver. 
In  our  results,  we  used  a  quasi-Newton  approach. 

Once  we  have  solved  for  js(t)  it  can  then  be  incorporated  into 
incremental  equations  for  film  resistance  and  capacity  loss.  It  is 
assumed  that js(t)  is  constant  over  some  small  time  interval  At,  and 
is  denoted  as  js[N]  for  the  Nth  interval.  Since,  according  to  Eq.(ll), 
the  change  in  film  thickness  is  proportional  to  js,  we  can  arrive  at 
an  incremental  equation  for  film  thickness  as 

4m[N]  =  4m [N  -  1]  -  ^Js[N  -  4 
anpt 

where  At  is  the  equation  update  period,  and  noting  that  the  sign  of 
js  is  negative.  This  result  can  be  used  to  calculate  the  film  resistance 
as 

RfilmlN]  =  RfilmlN  -  1  ]  -  -  1  ]. 

where  Requiv  =  W* Li  +  zLi2co3  /^u2co3  •  Similarly,  we  can  discretize 
the  capacity  equation  (13)  to  find  that 


QIN]  =  QIN  ~  1  ]  +  {ALnAtHslN  -  1  ]. 


3.6.  Summary  of  method 


4.  Simulation  and  results 


Here  we  summarize  the  complete  method  for  determining  the 
rate  of  lithium  loss.  The  method  is  iterative,  requiring  that  jn(t)  + 
js(t)  =  jtotal(t).  We  initialize  the  iteration  with  js(t)  =  0. 

1.  Compute  i0  n  based  on  the  present  cell  state  of  charge  and  set 

m  =  0. 

2.  Compute 


I(t) 

K&ALn  ’ 

3.  Set Ut)  =  (I(t)/(AL„))-Mt). 

4.  Compute 


E(t)  =  - 


/tec  -  (1  -  tp ) 

DfF 


+  1 


P (t)  =  +  ^^asinh 


m, 

6 


MO 

2anio,n 


+  Un  + 


jn  (QKfiln 
On 


The  validity  of  this  reduced-order  model  depends  first  on  the 
accuracy  of  the  underlying  partial  differential  equation  model, 
which  we  assume  here  to  be  exact.  It  then  depends  on  how 
closely  the  reduced-order  approximation  of  js  matches  the  exact 
calculation  of  js.  In  this  section,  results  from  both  the  full  and 
reduced-order  models  for  js  are  compared. 

To  compare  the  PDE  and  reduced-order  models,  we  conducted 
a  series  of  simulations.  In  each  simulation,  the  cell  was  initially  at 
rest.  A  sudden  pulse  of  current  was  then  applied,  and  the  result¬ 
ing^  from  the  PDE  model,  averaged  over  a  1-s  interval  subsequent 
to  the  pulse,  was  compared  to  the  computed  js  from  the  ROM.  To 
simulate  the  PDE  model,  we  used  COMSOL  Multiphysics  3.5a  [13] 
coupled  with  a  MATLAB[14]  script  to  cycle  through  the  series  of 
simulations  and  analyze  results.  Specifically,  each  simulation  com¬ 
prised  a  1.2  s  time  interval,  where  the  cell  current  iapp  was  modeled 
as  a  step  function,  which  was  applied  at  t  =  0.2  s.  We  found  that  the 
initial  rest  interval  facilitated  convergence  of  the  solution  by  allow¬ 
ing  the  PDE  solver  to  adjust  its  initial  conditions  before  applying  the 
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Table  1 

Electrode  parameters  for  simulation. 


Symbol 

Units 

Anode 

Separator 

Cathode 

L 

|jim 

85 

76.2 

179.3 

R 

|jim 

12.5 

- 

8.5 

A 

m2 

1 

1 

1 

a 

Sm_1 

100 

- 

3.8 

Ss 

- 

0.59 

- 

0.534 

£e 

- 

0.36 

1 

0.416 

Ke 

Sm_1 

0.2875 

0.2875 

0.2875 

brug 

- 

1.5 

- 

1.5 

rmax 

Ls 

mol  m3 

30  540 

- 

22  860 

c° 

mol  m3 

1000 

1000 

1000 

- 

0.10 

- 

0.95 

@i,  max 

- 

0.90 

- 

0.175 

Ds 

m2  s_1 

2.0  xlO-14 

- 

1.0  xlO”13 

De 

m2  s_1 

7.5  x  10-11 

7.5  x  10-11 

7.5  x  10-11 

t°+ 

- 

0.363 

0.363 

0.363 

k 

Am5/2  moP3/2 

2  x  10-6 

- 

2  x  10-6 

da 

- 

0.5 

- 

0.5 

dc 

- 

0.5 

- 

0.5 

da,s 

- 

0.3 

- 

- 

dc,s 

- 

0.7 

- 

- 

Uref,s 

V 

0.0 

- 

- 

Rsei 

Q,  m2 

0.002 

- 

- 

io,s 

A  ur2 

10 

- 

- 

step  current.  (Even  so,  convergence  of  the  PDE  simulations  proved 
troublesome,  and  required  user  vigilance  to  ensure  reliable  results.) 
The  solver  absolute  tolerance  was  set  to  10-1  and  the  relative  tol¬ 
erance  to  10-3.  The  default  direct  UMFPACK  solver  was  used,  with 
106  mesh  points  in  the  ID  battery  model  and  7056  elements  in 
the  2D  electrode  model.  The  solver  timestep  was  set  to  20  ms.  The 
ROM  results  were  computed  by  a  MATLAB  script  using  the  method 
summarized  in  Section  3.6. 


The  cell  parameters  that  we  used  in  the  simulations  match  those 
used  in  [9]  and  are  listed  in  Table  1 .  The  applied  current  was  varied 
between  0  C  and  3  C  in  increments  of  C/33;  the  initial  cell  SOC  was 
varied  between  0%  and  100%  in  steps  of  1%,  and  temperature  was 
held  constant  at  25  °C.  We  found  that  the  adjustable  tuning  factor 
/3  =  1.7  worked  well  (this  implies  the  change  in  electrolyte  concen¬ 
tration  near  the  separator  changes  nearly  twice  as  quickly  as  it  does 
near  the  current  collector).  A  total  of  10 100  simulations  were  run. 
As  one  point  of  comparison,  the  set  of  full-order  PDE  simulations 
took  approximately  12  h  to  complete,  utilizing  an  average  of  three 
cores,  on  an  Intel  i7  processor,  while  the  set  ROM  simulations  took 
approximately  21  s  to  complete,  utilizing  an  average  of  one  core 
on  the  same  machine.  The  speedup,  on  a  per-simulation  per-core 
basis,  is  more  than  5000 : 1.  This  is  the  primary  advantage  of  the 
ROM  over  the  PDE  model. 

Lithium  plating  occurs  when  the  side-reaction  overpotential  is 
negative  {rjs  <  0).  Fig.  2(a)  illustrates  the  side-reaction  overpotential 
across  the  anode  for  this  cell  model,  where  x  =  0  is  adjacent  to  the 
current-collector  and  x  =  85  [xm  is  adjacent  to  the  separator,  imme¬ 
diately  following  the  onset  of  a  charge  current  pulse.  From  the  PDE 
result,  we  expect  lithium  deposit  to  occur  between  about  x  =  42  |xm 
and  x  =  85  pan.  From  the  ROM,  we  expect  lithium  deposit  to  occur 
between  x0  =  49[xm  and  x  =  85[xm.  Fig.  2(b)  shows  the  resulting 
rate  of  lithium  deposition  for  the  PDE  and  ROM  solutions.  The  time- 
average  deposition  rate  of  the  ROM  is  somewhat  higher  than  the 
time-average  deposition  rate  of  the  PDE  over  the  1  s  interval. 

Fig.  3  illustrates  the  predicted  overcharge  rates  over  all  scenar¬ 
ios  for  the  PDE  and  the  ROM  solutions.  As  expected,  deposition  is 
worse  at  high  SOC  and  high  charge  rates.  The  PDE  and  ROM  solu¬ 
tions  generally  agree  very  well,  with  greatest  mismatch  at  high 
charge  rates. 

Fig.  4  shows  a  different  view  of  Fig.  3.  Cross  sections  through 
both  the  PDE  and  ROM  solution  spaces  are  plotted  and  compared. 


(a)  Overpotential  predictions  at  90%  SOC,  2.0C  rate 


(b)  Overcharge  predictions  at  90%  SOC,  2.0C  rate 


Fig.  2.  Comparing  PDE  and  ROM  overpotential  and  deposition-rate  predictions. 


Fig.  3.  Comparing  the  PDE  and  ROM  over  multiple  scenarios. 
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(a)  Instantaneous  degradation  rate  as  function  of  SOC  (b)  Instantaneous  degradation  rate  as  function  of  current 


Fig.  4.  Cross  sections  through  the  plots  in  Fig.  3  to  more  clearly  illustrate  effects  of  SOC  and  rate. 


(a)  Regions  where  models  predict  overcharge 


Modeling  error  between  PDE  and  ROM  (mA  cm'3) 


Charge  current 


Fig.  5.  Comparing  the  PDE  and  ROM  solutions. 


Frame  (a)  shows  how  the  two  methods  compare  where  each  pair  of 
lines  represents  a  specific  charge  rate.  As  noted  before,  but  perhaps 
more  clearly  seen  here,  the  difference  between  the  PDE  and  ROM 
solutions  are  greatest  at  high  charge  rates.  Frame  (b)  shows  how  the 
two  methods  compare  where  each  pair  of  lines  represents  a  specific 
initial  SOC.  The  difference  is  greatest  at  moderate  SOC  levels. 

Finally,  Fig.  5  illustrates  the  error  between  the  PDE  and  ROM 
solutions  in  two  ways.  Frame  (a)  shows  the  regions  where  the  two 
methods  agree  on  whether  lithium  deposition  will  occur,  and  the 
region  where  they  disagree.  The  region  of  disagreement  is  the  very 
narrow  sliver  at  around  2.4  C  and  25%  SOC,  where  the  ROM  pre¬ 
dicts  overcharge  but  the  PDE  does  not.  Otherwise,  the  boundaries 
are  identical.  Frame  (b)  shows  the  error  between  the  solutions, 
calculated  as  js  R0M  -js, pde-  The  maximum  error  is  approximately 
65mA  cm-2.  When  compared  to  Fig.  3,  we  see  that  the  relative  error 
is  on  the  order  of  10%. 

For  the  purpose  of  control  system  design,  the  results  of  Fig.  5(a) 
are  the  most  important.  Since  lithium  deposition  is  such  a  severe 
degradation  mechanism,  a  charging  control  scheme  should  avoid 
ever  commanding  a  control  action  that  would  cause  any  lithium 
deposition  to  occur.  A  time-optimal  charger,  based  only  on  the  PDE 
model  of  lithium  deposition,  would  select  charge  pulse  current  to 
follow  the  upper  contour  in  Fig.  5(a).  This  allows  the  maximum 
charge  rate  at  any  point  in  time,  while  causing  no  lithium  plat¬ 
ing.  In  comparison,  a  time-optimal  charger,  based  only  on  the  ROM 
model  of  lithium  deposition,  would  select  charge  pulse  current  to 
follow  the  lower  contour  in  the  figure.  This  will  result  in  somewhat 
slower  charging.  But,  because  the  ROM  over-predicts  the  amount 
of  lithium  deposition,  it  will  also  result  in  a  charging  scheme  that 
is  conservative,  which  is  a  beneficial  feature. 

We  conducted  additional  simulations  to  investigate  the  effect 
of  pulse  length  At  in  Assumption  1  of  Section  3.  That  is,  how  long 
can  the  charge  pulse  be  before  the  full-order  PDE  model  and  the 


reduced  order  model  results  are  significantly  different?  We  found 
that  pulse  lengths  less  than  10  s  are  generally  well  matched,  but 
pulse  lengths  much  greater  than  10s  can  give  significant  PDE  versus 
ROM  mismatch.  For  long  pulse  durations,  the  quasi-static  nature 
of  Assumption  1  is  violated,  and  an  a  significant  offset  is  noted  in 
actual  time-varying  0e  versus  the  at-rest  0e,  moving  the  crossover 
point  of  r]s{x,  t).  This  causes  the  ROM  to  under  predict  the  value 
of  lithium  plating  computed  by  the  PDE.  For  this  reason,  we  pro¬ 
pose  that  the  ROM  is  of  most  value  for  computing  current  limits  in 
dynamic  applications  such  as  hybrid-electric  vehicles,  where  a  bias 
in  4>e  cannot  develop  due  to  the  random  nature  of  power  demand, 
but  is  of  less  value  for  controlling  full  charges,  such  as  for  electric 
vehicle  applications.  A  future  paper  will  disclose  more  advanced 
reduced-order  methods  to  compute  lithium  plating  for  prolonged 
charge  events  as  well  as  short  charge  events. 

We  make  one  final  comment  regarding  efficiency.  The  speedup 
of  the  ROM  versus  PDE  solutions  can  be  much  greater  than  the 
aforementioned  5000:1  if  the  ROM  solutions  are  pre-computed 
and  stored  in  a  table  and  accessed  via  table  lookup.  Subsequently, 
“computing”  any  value  of  js  R0M  would  be  nearly  instantaneous.  We 
note  that  the  ROM  solution  changes  as  the  film  resistance  changes, 
but  the  film  resistance  changes  very  slowly.  The  entire  table  might 
be  updated  by  the  battery  management  system  once  per  opera¬ 
tional  period  (e.g.,  once  per  day),  and  then  utilized  throughout  that 
operational  period  for  significant  performance  gains. 

5.  Conclusions 

From  the  results  above,  we  see  excellent  correlation  for  the  rate 
of  lithium  deposition  on  overcharge  between  a  full-order  coupled 
multi-physics  PDE  model  and  the  reduced-order  model  derived  in 
this  paper.  The  reduced  order  model  takes  the  very  difficult  and 
time  consuming  procedure  to  calculate  the  governing  equations 
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necessary  to  obtain  js  in  a  PDE  model,  and  reduces  them  to  a  very 
fast  algebraic  procedure  that  a  small  microcontroller  can  easily  cal¬ 
culate.  As  noted,  the  speedup  is  at  least  5000:1  when  using  the 
ROM  versus  the  PDE  solver,  and  can  be  made  much  faster  via  the 
proposed  table-lookup  approach  to  storing  results. 

Of  particular  importance,  the  reduced-order  model  is  able  to 
closely  but  conservatively  predict  the  boundary  between  when 
lithium  plating  will  not  and  will  occur  in  the  SOC/charge-current 
space,  especially  for  short-duration  charge  events.  This  boundary 
is  critical  for  producing  controls  that  will  not  prematurely  age  a 
cell  due  to  overcharge.  In  a  future  paper,  the  authors  will  explain 
their  current  research  into  optimal  charging  profiles  that  consider 
overcharge  as  well  as  other  degradation  mechanisms. 

Improvements  to  the  ROM  might  be  possible  if  the  assumptions 
made  in  Section  3  can  be  relaxed.  In  particular,  Assumption  3  is  not 
very  accurate  at  high  rates  of  charge.  However,  we  find  the  present 
results  to  be  adequate  for  initial  development  of  control  algorithms 
for  battery  charging. 

List  of  symbols 


a  specific  surface  area  of  porous  electrode  (m-1 ) 

c  concentration  of  Li  or  Li+  ions  (mol  m-3 ) 

D  diffusion  coefficient  (m2  s-1 ) 

E  composite  quantity  defined  to  ease  notation  per  Eq.  (19) 

(Vm-2) 

F  Faraday’s  constant  (96  487  C  mol-1 ) 

ie  current  density  of  electrolyte  (A  nrr 2 ) 

i0  exchange-current  density  for  intercalation  reaction 

(Am-2) 

i'se  current  density  of  solid-electrolyte  interface  (A  m-2 ) 
i0,s  exchange-current  density  for  side  reaction  (A  m-2 ) 
j  local  volumetric  current  density  for  intercalation  reaction 

(Am-3) 

js  local  volumetric  current  density  for  side  reaction  (A  m-3 ) 

k  rate  constant  of  electrochemical  reaction  (A  m5/2  mol-3/2 ) 

L  length  of  cell  or  electrode  (m) 

P  composite  quantity  defined  to  ease  notation  per  Eq.  (23) 

(V) 

Q.  cell  capacity  (C) 

R  particle  radius  (pan) 

Rfilm  film  resistance  at  the  electrode/electrolyte  interface 

(£2m2) 

Rg  universal  gas  constant  (8.314  J  mol-1  K-1 ) 

T  temperature  (K) 

t  time  (s) 

transport  number 

U  local  equilibrium  potential  (V) 

x  length  dimension  (m) 


Greek 

ota,  otc  anodic  and  cathodic  transfer  coefficients  of  electrochem¬ 
ical  reaction 

e  volume  fraction  of  a  phase 

0  local  potential  of  a  phase  (V) 

1 7  local  overpotential  driving  electrochemical  reaction  (V) 

k  conductivity  of  electrolyte  (S  m-1 ) 

0  state-of-charge  or  stoichiometry  of  electrode 

a  conductivity  of  electrode  (S  m-1 ) 

p  density  of  active  material  (kg  m-3 ) 

8  film  thickness  (m) 

Subscript/superscript 

e  pertaining  to  the  electrolyte  phase 

s  pertaining  to  the  solid  phase 

n  pertaining  to  the  negative  electrode 

p  pertaining  to  the  positive  electrode 
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